Parameter-free dissipation in simulated sliding friction 
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Non-equilibrium molecular dynamics simulations, of crucial importance in sliding friction, are 
hampered by arbitrariness and uncertainties in the way Joule heat is removed. We implement in 
a realistic frictional simulation a parameter-free, non-markovian, stochastic dynamics, which, as 
expected from theory, absorbs Joule heat precisely as a semi-infinite harmonic substrate would. 
Simulating stick-slip friction of a slider over a 2D Lennard- Jones solid, we compare our virtually ex- 
act frictional results with approximate ones from commonly adopted empirical dissipation schemes. 
While the latter are generally in serious error, we show that the exact results can be closely repro- 
duced by a viscous Langevin dissipation at the boundary layer, once the back-reflected frictional 
energy is variationally optimized. 

PACS numbers: 81.40.Pq, 68.35. Af, 05.10.Gg, 46.55.+d 
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The role of Molecular Dynamics (MD) simulation in 
the theory of sliding friction and nanofriction can hardly 
be overestimated [l|, [2[ . In a typical simulation, a slider 
is driven by an external force over a simulated solid sub- 
strate whose atoms, interacting through realistically cho- 
sen interatomic forces, vibrate and move according to 
Newton's law. Naturally, in order to attain a frictional 
steady state the Joule heat must be removed. Unfortu- 
nately, a realistic energy dissipation is generally impos- 
sible to simulate reliably, owing to size (and time) limi- 
tations. In a small size simulation cell, the phonons gen- 
erated at the sliding interface are back-reflected by the 
cell boundaries, rather than propagated away to prop- 
erly disperse the Joule heat. The ensuing problem is a 
spuriously accumulating phonon population in the slider- 
substrate interface region. The empirical introduction 
in the equations of motion of ad-hoc Langevin viscous 
damping terms —mjiqi (with m and qi the mass and the 
velocity of the i-th particle) and of an associated random 
noise, corresponding to some "thermostat" temperature 
T represents the handiest and commonest solution. 
However, both this procedure and the choice of ther- 
mostat and damping parameters ji are vastly arbitrary. 
The problem is not just one of principle. Dry friction 
generally involves stick-slip [l[, and each slip generates 
a burst of phonons of very specific nature and compo- 
sition. In realistic simulations the overall effect on fric- 
tional dynamics of the partial back-reflection of this burst 
must be minimized. Unless dealt with, back-reflection 
will cause the simulated steady state and friction coeffi- 
cient to depend upon unphysical damping parameters. In 
the Prandtl-Tomlinson model 0] for instance, damping 
is known to modify the kinetics and the friction, both in 
the stick-slip regime (including the multiple-slips seen in 
Atomic Force Microscopy [f|) and in the smooth sliding 
state. To a lesser or larger extent this lamentable state 
of affairs is common to all MD frictional simulations. 

Pursuing a viable solution, one wishes to modify the 
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FIG. 1: (Color online) (a) Sketch of the simulated sys- 
tem. (b),(c) Layer-averaged kinetic energy time evolution 
of a surface-injected phonon burst. Phonon back-reflection 
is strong without dissipation at the boundary layer (c), but 
accurately canceled once the correct dissipative kernels are 
included (b). 



equations of motion inside a relatively small simulation 
cell, so that they will reproduce the frictional dynamics 
of a much larger system, once the remaining variables are 
integrated out. Integrating out degrees of freedom is a 
classic problem, largely analysed in literature d,|B-[i|. In 
the context of MD simulation, Green's function methods 
were formulated for quasistatic mechanical contacts [l(| ; 
approaches based on a discrete-continuum matching have 



also been discussed [ll|, |12(. Among others, a formally 



exact dissipation formulation was given early on by Ru- 
bin Q . Integrating out N — 1 atoms in a linear harmonic 
chain yields a non-markovian Langevin-like equation of 
motion of a single atom of interest. Extensions of this 
method were applied to a variety of problems, including 
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the relaxation of impurity molecules in solids 13j, Il4j . 
and atom-surface scattering flBl . [l6| . Detailed formula- 
tions suitable for MD simulations were given by Li et 
al. and by Kantorovich 0, H Qj}- The idea to vari- 
ationally minimize back-reflection of individual phonons 
has also been put forward 0. However, this body of the- 
ory has not yet found its way into realistic sliding friction 
MD simulations, so that neither the importance of back- 
reflection errors in dry friction, nor a realistic way to get 
rid of them have really been demonstrated. Here we de- 
scribe how both goals are achieved, in a simulated 2D 
Lennard- Jones system that exhibits a realistic stick-slip. 
We find that exact stick-slip friction is indeed different 
from that of empirical damping schemes. However the 
back-reflected energy can be variationally minimized - 
although differently from Li et al. — the resulting ap- 
proximate dynamics now reproducing quite accurately 
the exact benchmark stick-slip friction. 

We focus on a model sliding system made of a semi- 
infinite two-dimensional crystalline substrate with N z 
monatomic layers each of N x atoms, and of a slider made 
of a single chain of N' x atoms. All atoms interact via 
a Lennard- Jones (LJ) potential (cutoff for simplicity to 
first- neighbors). The slider is pressed against the sub- 
strate by a normal "load" force F 0: and is driven along 
x (parallel to the substrate) through a spring of strength 
fc, whose end is pulled at constant velocity v . Following 
similar earlier formulations 0-0, the ideal infinitely thick 
substrate is divided as in Fig.QJa), into three regions: (i) 
a live slab comprising N z — 1 atomic layers whose motion 
is fully simulated by Newton's equations; (ii) the dissi- 
pative boundary (A^-th) layer, whose motion includes 
the effective non-markovian Langevin terms; (iii) the re- 
maining semi-infinite solid acting as a phonon absorber, 
or heat bath, whose degrees of freedom are integrated 
out, providing the source of effective damping terms in 

(ii) . In our practical implementation we substitute the 
full LJ potential within region (iii) and between (ii) and 

(iii) with its harmonic approximation. This approxima- 
tion, necessary in order to get an analytical form for the 
effective forces on the boundary atoms, is all the more ac- 
curate the weaker the intensity of frictional phonons. In 
practice, for crystal substrates above quantum freezing 
and below the Debye temperature, this level of accuracy 
can always be attained by a sufficient thickness N z — I 
of the live simulation cell (i). The harmonic heat bath 
(iii) is decoupled by diagonalizing the dynamical matrix 
D^ v (k,l denoting atoms), obtaining eigenvalues and 
eigenvectors A^. The equations of motion for the bound 
ary layer atoms are 0, 
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where i and j denote boundary layer atoms, \i and 
v indicate x/z components, U su b is the LJ interaction 
between the i-th boundary atom and those inside the 
simulation cell (i). The second term is non-markovian 
and non-conservative, introducing an effective damping 
proportional to the velocity of the j-th atom, through 
a time convolution with the memory kernel functions 
K l j v (t). Standard kernels are built from the harmonic 
eigenvalues and eigenvectors of the heat-bath dynami- 
cal matrix and from the coupling vectors </>^ containing 
the harmonic coupling constants £)?* and D l * v of the 
i-th atom of region (ii) with the fc-th heat-bath atom 

(A»<X*i-A„ cos(wfet) . They oscil- 



late and decay with time, with power law tails due to 
the bath acoustical phonon branches. Periodic boundary 
conditions along the x direction guarantee translational 
invariance, so that K l ^ v {t) is a function of \i — j\ only. 
As kernels inherit their symmetry properties from those 
of the heat-bath dynamical matrix, one can show that 
K%{t) = -Kg(t) and K%{t) = K^(t). As \i - j\ 
grows, \KJjj v (t)\ decrease, but again not exponentially, 
so that correlations must be included up to large dis- 
tance 7] . By cutting kernels off at time r c one could limit 
the time-integrals in Eq. ([1} , which need to be calculated 
at each time step. By increasing both r c and the live 
simulation cell size N z , one refines the dissipation scheme 
accuracy as much as desired, although at increased com- 
putational cost. The third term in Eq. (TTJ) is a gaus- 
sian stochastic force present at non-zero temperature, re- 
sponsible for the energy transfer between the heat-bath 
(iii) and the live substrate (i), with {R l ^(t)) = 0, and 
(Ri(t)Ri(t')) = mk B TK l ^(t - t') where brackets de- 
note ensemble average, ks is Boltzmann's constant and 
T the temperature. As is well known 0], this relationship 
fulfils the fluctuation-dissipation theorem — the scheme 
represents a well defined thermostat in contact with the 
simulation cell. The last term in Eq. (fTJ) is the harmonic 
coupling between atoms i-th and j-th within the bound- 
ary layer, where the coupling constant D % j v is modified 

byi^i(o). 



As a first step before simulating friction we implement 
the set of equations |T]), along with ordinary Newton's 
equations governing the live cell in an MD simulation, 
to test phonon reflection. Starting with the equilibrated 
system (N z = 30 close packed layers and N x — 10 atoms 
per layer) at temperature T, a test burst of phonons is 
introduced at the topmost layer and time-evolved. Re- 
sults in Fig. [IJb) show how a relatively thin N z = 30 
layer substrate (i+ii) is able to mimick, within the ex- 
act dissipation scheme, the full ideal semi-infinite system 
(i+ii+iii). Layer-resolved kinetic energies inside the sim- 
ulated substrate show the group of phonons propagating 
below the surface. Upon reaching the boundary layer 
the phonons are perfectly absorbed as if they propagated 




FIG. 2: (Color online) (a) Simulated friction force F(t) for 
the "exact" non-markovian dissipation scheme of Eq. fl]), and 
[(b), (c) and (d)] for empirical Langevin schemes with a —jq 
damping applied to tip LI, whole substrate L2 and bottom 
layer L3. The 7 value is identified by numbers 1-5 in Fig[3] 
Dashed lines: mean value (F). 



into the (integrated out) semi- infinite crystal (hi). For 
comparison, Fig. QJc) shows the same phonons massively 
back reflected once the memory kernels are removed from 
the boundary layer. 

We next simulate stick-slip sliding friction by driving 
a slider, here consisting of a LJ chain of N' x = 9 atoms, 
over the same live substrate as above |18| . The instanta- 
neous friction force is measured by the spring elongation 
F(t) = k(xcM — Vot), xcm being the slider center of mass 
position. The sawtooth force profile typical of stick-slip 
friction is obtained (Fig. [2] (a)). The friction coefficient, 
obtained by averaging over a stationary sequence of stick- 
slips is (F)/F = 0.116 ± 0.002. The stick-slip pattern 
is irregular, with a periodicity similar but not exactly 
matching, a substrate lattice spacing. The highest spikes 
signal the forward jump of most slider atoms, smaller 
ones involve only about 2/3 of them. A measure of the 
distribution of the spike heights is the variance of F(t), 
i.e., a — - } F \. 2 J^ s [F(t) — (F)] 2 dt, where r s is the total 



simulation time. These simulations using the full Eq. (fl). 
and the corresponding frictional results can be considered 
essentially exact, and represent our benchmark reference 
of stick-slip friction with a correct Joule heat removal. 
The numerical implementation of this standard scheme 
in a generic 3D case will in principle be possible but 
time consuming, and far from handy. Long-range non- 
markovian correlations imply a general scaling as N% , so 
that simulations for large-size 3D sliding systems may 
pose a practical challenge of parallel computing. 

The virtually exact frictional dynamics just obtained 
can now be compared with empirical Langevin schemes 
commonly adopted in friction simulations, where a vis- 
cous damping —75 is arbitrarily applied to the motion 
of some atoms in the system, for instance to slider 
atoms (LI, dashed line in FigOJ) or to all substrate 
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FIG. 3: (Color online) (a),(b) Friction coefficient (F)/F and 
variance (a) as a function of a viscous damping 7 for different 
empirical Langevin dissipation schemes, compared with ex- 
act non-markovian values (blue stripes) . (c) Average internal 
energy W per substrate atom for the markovian thermostat 
L3 divided by the internal energy Wexact per substrate atom 
for the non-markovian scheme {TJ. The exact and empirical 
frictional behaviour nearly coincide when the boundary layer 
parameter 7 is such as to variationally minimize the relative 
total back-reflected energy W/W exa ct- 



atoms (L2, dotted line), or to the boundary atomic 
layer only (L3, solid line). As anticipated, the stick- 
slip friction simulated with all these empirical schemes 
depends vastly on 7, and generally deviates systemat- 
ically from the correct benchmark. The closest agree- 
ment is obtained in case L3 where the viscous damping 
—mjq^t) is applied to the i-th boundary layer atom mo- 
tion, mql(t) = - rwyqlit) + f2*(i), with the appro- 



priate gaussian stochastic force Ri(t) with (£)) = 
and (Ri(t)Ri(t')) = 2mk B T 1 d^8 lJ d(t - t'). Since 
the parameter 7 is adjustable, we may seek to opti- 
mize it in order to approximate the exact results. Li 
et at 7] considered variationally minimizing a group- 
velocity weighted phonon reflection. In the actual fric- 
tion simulation, the steady state internal energy increase 
W = (E) — (Eq) over the equilibrium value Eq is eas- 
ily computable. Our proposed scheme is to variationally 
minimize W, a quantity generally positive and propor- 



4 



0.12 



0.11 



0.10 



0.09 




20 30 



FIG. 4: (Color online) Friction coefficient as a function of 
the simulated cell thickness N z for different driving velocities 
vo = 0.01 (blue) and vo = 0.02 (black) obtained with the 
non-markovian dissipation scheme. The inset reports the av- 
erage friction force for different loads Fo obtained with the 
non-markovian approach (blue stripes) and the comparison 
with the markovian Langevin scheme L3 at different 7 (black 
curves) . 



tional to the Joule frictional heat. In absence of back- 
reflection W is minimal, whereas partial energy reflec- 
tion at the boundary will cause Joule heat to artificially 
accumulate in the slab and increase W, so that for any ar- 
bitrary dissipation scheme W/W exac t > 1- As a function 
of 7, a variational minimum of W/W exac t occurs because 
back-reflection of frictional phonons is large both when 
the boundary layer damping 7 is too small and too large. 
This ratio is shown in Fig. OHc). The agreement between 
the exact frictional results, where no phonons are back 
reflected, and the variational one at 7 = 10 is excellent. 
(We note incidentally that at the optimal 7 the friction 
coefficient also peaks). To confirm the variational result, 
we changed system parameters, including sliding veloc- 
ity, and load. The variable load results of Fig. 0] (inset) 
show that the coincidence of optimal and exact friction 
is systematic. Fig. 2] also shows a dependence of calcu- 
lated friction upon the thickness of the simulated sub- 
strate portion (i+ii), converging for sufficiently large N z . 
The required thickness of a few tens of layers depends on 
velocity, and is reached once the substrate inertia grows 
large enough to stop interfering spuriously with the fric- 
tional dynamics. 

In conclusion, we demonstrated frictional MD sim- 
ulations with a nonmarkovian dissipation, enacting 
the correct disposal of generated phonons, even in the 
relatively violent stick-slip case. Using the virtually 
exact frictional results so obtained as a reference, we 
benchmarked common empirical viscous dissipation 
schemes, finding them generally wanting. There exists 
however an optimal viscous damping j opt which, once 



applied to the cell boundary layer, yields results nearly 
indistinguishable from the exact ones. The optimal 
damping parameter variationally minimizes the total 
back-reflected energy at the cell boundary, and can be 
identified even without any exact reference calculation. 
This optimal damping scheme, which unlike the exact 
one does not require a knowledge of the substrate 
vibrational properties, is a good candidate for adoption 
in future practical MD frictional simulations. 
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Sliding simulations are performed at temperature fcsT = 
0.035, roughly corresponding to T/T meit i n9 = 0.06 (LJ 
units used throughout). The load is Fo — 10, the av- 
erage sliding velocity vo = 0.01, and the spring con- 
stant k = 5. Periodic boundary conditions are applied 
to substrate and slider along x. To favor sliding, the 
strength of the slider-substrate LJ interaction is reduced 
from 1 to 0.6. Equations of motion are integrated by a 
modified velocity- Verlet algorithm with a time step of 
At = 5 • 10~ 3 . Memory kernels are cutoff at r c = 5 • 10 3 



time-steps, with good accuracy of the frictional pattern 
in all its aspects, including phonon dissipation, our main 



